# -*- coding: utf-8 -*-

"""
Created on Sun Feb  4 18:39:12 2024

@author: Nicholas Drachman
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['svg.fonttype'] = 'none'

import matplotlib.gridspec as gridspec



''' Fig 4b,c'''

pathLVM = 'Fig 4/b/exp 4.lvm'
df = pd.read_csv(pathLVM, header = 21, delimiter = '\t')

Vt = df['Tip voltage']
Ie = df['Emission current']
It = df['Detected current']
Eff = df['Efficiency']

tmax = len(df)/2
t = np.linspace(0,tmax,len(df))


tRange = np.array([3650,5800])

Ie_new = Ie - 1e-10  # Measured offset in emission current (attributed to leakage through coax)
Eff_new = It/Ie_new


fig = plt.subplots()
# Define the grid
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])  # 2x1 grid, top plot is twice the height of the bottom plot
# Create the subplots
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])

ax0.plot((t[tRange[0]:tRange[1]] - t[tRange[0]])/60, Ie_new[tRange[0]:tRange[1]]*1e12, linewidth = 1, label = 'Emitted current')
ax0.plot((t[tRange[0]:tRange[1]] - t[tRange[0]])/60, It[tRange[0]:tRange[1]]*1e12, linewidth = 1, label = 'Transmitted current')
ax0.set_xlim(0,14)
ax0.set_ylabel('Current (pA)', fontsize = 8, labelpad = 0)
ax0.tick_params(which = 'both', direction = 'in', labelsize = 8)
ax0.set_ylim(0, 1200)
ax0.legend()

ax1.plot((t[tRange[0]:tRange[1]] - t[tRange[0]])/60, 100*Eff_new[tRange[0]:tRange[1]], 'k')
ax1.set_ylim(0,100)
ax1.set_xlim(0,14)
ax1.set_ylabel('efficiency (%)', fontsize = 8, labelpad = 1)
ax1.set_xlabel('time (min)', fontsize = 8, labelpad = 1)
ax1.tick_params(labelsize = 8, direction = 'in')


plt.tight_layout(w_pad = 0.05, h_pad =0.3, pad = 0)





''' Fig 4d,e '''

pathLVM = 'Fig 4/c/Exp 16.lvm'

# Read the data
df = pd.read_csv(pathLVM, header = 21, delimiter = '\t')
df.head()
Vt = df['Tip voltage']
Ve = df['Extractor Voltage']
Ie = df['Emitted Current']
It = df['Tramsmitted current']
Eff = df['Efficiency']
t = np.array(df['Time'])

t = (t-t[0])/1e3


import matplotlib.gridspec as gridspec
from scipy.signal import bessel
from scipy.signal import filtfilt

a,b = 4375,4860   # timerange to plot


Ie_offset = 4.9e-12  # Measured offset in emission current measurement
It_offset = 2e-13   # Measured offset in transmitted current measurement

Ie_data = (Ie[a:b] - Ie_offset)*1e12
It_data = (It[a:b] - It_offset)*1e12

time = t[a:b] - t[a]

sampling_interval = np.median(np.diff(time))
sampling_frequency = 1 / sampling_interval
print(sampling_frequency)

# Design a Butterworth low-pass filter
cutoff_frequency = 0.1  # Cutoff frequency in Hz
order = 8 # Filter order
bb, aa = bessel(order, cutoff_frequency / (0.5 * sampling_frequency), btype='low')

# Apply the filter
Ie_filt = filtfilt(bb, aa, Ie_data)
It_filt = filtfilt(bb, aa, It_data)


Eff_filt = It_filt/Ie_filt

fig = plt.figure()

# Define the grid
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])  # 2x1 grid, top plot is twice the height of the bottom plot
# Create the subplots
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])


ax0.plot(time / 60, Ie_filt, lw = 0.75, label = 'Emitted current')
ax0.plot(time / 60, It_filt, lw = 0.75, label = 'Trasmitted current')
ax0.set_ylim(0,8)
ax0.set_xlim(0,14)
ax0.legend()
ax0.set_ylabel('current (pA)', fontsize = 8, labelpad =5)
ax0.tick_params(labelsize = 8, direction = 'in')



ax1.plot(time/60, Eff_filt*100, color = 'k', lw = 1, label = 'Efficiency')
ax1.set_ylim(0,100)
ax1.set_xlim(0,14)
ax1.set_ylabel('efficiency (%)', fontsize = 8, labelpad = 0.5)
ax1.set_xlabel('time (min)', fontsize = 8, labelpad = 1)
ax1.tick_params(labelsize = 8, direction = 'in')
plt.tight_layout(w_pad = 0.05, h_pad =0.15, pad = 0)




''' Fig 4g,h '''

path = 'Fig 4/e/exp2.lvm'

df = pd.read_csv(path, header = 21, delimiter = '\t')
    
df.head()

Vt = df['Tip voltage']
Ie = df['Emission current']
Id = df['droplet current']
Ii = df['ion current']

tmax = len(df)/2
t = np.linspace(0,tmax,len(df))

ion_frac = Ii/(Ii+Id)

trange = np.array([2260,2388])*2

fig = plt.subplots()

# Define the grid
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])  # 2x1 grid, top plot is twice the height of the bottom plot
# Create the subplots
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])

ax0.plot()


ax0.plot((t[trange[0]:trange[1]] - t[trange[0]])/60,Ii[trange[0]:trange[1]]*1e12, label = 'Ion current')
ax0.plot((t[trange[0]:trange[1]] - t[trange[0]])/60,Id[trange[0]:trange[1]]*1e12, label = 'Droplet current')
ax0.set_ylabel('current (pA)', fontsize = 8, labelpad = -0.5)
ax0.tick_params(labelsize = 8, direction = 'in')
ax0.set_ylim(-1,125)
ax0.set_xlim(0,2.1)
ax0.legend()

ax1.plot((t[trange[0]:trange[1]] - t[trange[0]])/60,(Ii/(Ii+Id))[trange[0]:trange[1]]*100, color = 'k')
ax1.set_xlim(0,2.1)
ax1.set_ylim(0,105)
ax1.set_xlabel('time (min)', fontsize = 8)
ax1.set_ylabel('ion fraction (%)', fontsize = 8, labelpad = 0)
ax1.tick_params(labelsize = 8, direction = 'in')
ax1.tick_params(axis = 'y', pad = 2)

plt.tight_layout(w_pad = 0.05, h_pad =0.3, pad = 0)


''' Fig 4i,j '''


path = 'Fig 4/f/Exp 5.lvm'

df = pd.read_csv(path, header = 21, delimiter = '\t')
    
Vt = df['Tip voltage']
Ie = df['Emission Current']
Ii = df['Ion Current']
Id = df['Droplet Current']
Eff = df['Efficiency']
Ve = df['Extractor Voltage']

t = np.array(df['Time'])

t = (t-t[0])/1e3


a,b = 3020,3300

N = 1
Ie_smooth = np.convolve(Ie,np.ones(N)/N, mode = 'valid') 
Ii_smooth = np.convolve(Ii,np.ones(N)/N, mode = 'valid')
Id_smooth = np.convolve(Id,np.ones(N)/N, mode = 'valid')
t_smooth = np.convolve(t,np.ones(N)/N, mode = 'valid')

Ie_offset = 4e-12
Ii_offset = .3e-12
Id_offset = 1.3e-13

Idet = Ii_smooth + Id_smooth - Ii_offset - Id_offset  # Total detected current

Eff_smooth = Idet/(Ie_smooth - Ie_offset)

ion_frac = (Ii_smooth - Ii_offset) / (Idet)



fig = plt.subplots()

# Define the grid
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])  # 2x1 grid, top plot is twice the height of the bottom plot
# Create the subplots
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])


# ax0.plot((t_smooth[a:b]-t_smooth[a]) / 60, (Ie_smooth[a:b] - Ie_offset)*1e12, label = 'emitted current')
ax0.plot((t_smooth[a:b]-t_smooth[a]) / 60, (Ii_smooth[a:b] - Ii_offset)*1e12, lw = 1, label = 'Ion current')
ax0.plot((t_smooth[a:b]-t_smooth[a]) / 60, (Id_smooth[a:b] - Id_offset)*1e12, lw = 1, label = 'Droplet current')
ax0.set_ylim(-0.02,3) 
ax0.set_xlim(0, 7.5)
ax0.set_ylabel('current (pA)', fontsize = 8, labelpad = 5)
ax0.legend()
ax0.tick_params(labelsize = 8, direction = 'in')


ax1.plot((t_smooth[a:b]-t_smooth[a]) / 60,ion_frac[a:b]*100, 'k')
ax1.set_xlim(0,7.5)
ax1.set_ylim(0,105)
ax1.set_xlabel('time (min)', fontsize = 8)
ax1.set_ylabel('ion fraction (%)', fontsize = 8, labelpad = 0)
ax1.tick_params(labelsize = 8, direction = 'in')
ax1.tick_params(axis = 'y', pad = 2)
plt.tight_layout(w_pad = 0.05, h_pad =0.3, pad = 0)

